VARIABLE ELIMINATION IN POST-TRANSLATIONAL 
MODIFICATION REACTION NETWORKS WITH MASS-ACTION 

KINETICS 



^^ Abstract. We define a subclass of Cfiemical Reaction Networks called Post-Translational 

f^ Modification systems. Important biological examples of such systems include MAPK 

cascades and two-component systems which are well-studied experimentally as well as 
theoretically. The steady states of such a system are solutions to a system of polynomial 
equations with as many variables as equations. Even for small systems the task of finding 
the solutions is daunting. We develop a mathematical framework based on the notion of 
a cut, which provides a linear elimination procedure to reduce the number of variables 
in the system. The steady states are parameterized algebraically by a set of "core" vari- 
ables, and the non-negative steady states correspond to non-negative values of the core 
variables. Further, minimal cuts are the connected components in the species graph and 
provide conservation laws. A criterion for when a set of independent conservation laws 
can be derived from cuts is given. 
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1. Introduction 



Signaling systems play an important role in regulation of cellular processes and are essen- 
tial for cellular decision making. Typical signaling systems react to stimulus in the (cellular) 
environment and transmit a signal through connected layers of biochemical species. The 
C<^ layers provide means to adjust the response according to the stimulus. A common form 

[^^ of signaling systems is Post-translational Modification (PTM) systems where species are 

^> activated in chemical reactions in order to propagate the signal through the system. 

'""' PTM systems have attracted considerable theoretical attention due to their abundance 

• • in nature [12] and regular form [17j . The dynamics can be modeled as -^^ = p{x) , where 

.!^ X = {xi, . . . ,Xn) are the variables (concentrations of species) of the system and p{x) is 

^ a vector of polynomials in x. Only certain types of reactions are allowed, restricting the 

form of p(x). In particular, small specific systems have been scrutinized, focusing on the 
dynamical behavior and the steady states of the systems. Examples include the biologically 
important MAPK cascades [121 [131 US] , as well as simpler signaling cascades [9l [TTl [20] . 

We focus on the steady states of a PTM system (defined formally in the next section) 
and how to determine them. Taken with mass-action kinetics, the system's steady states 
are solutions to a set of polynomial equations in the species and with coefficients given by 
unknown kinetic rates (i.e. unspecified parameters). In particular, the number of equations 
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to be solved is equal to the number of species. Even small systems might have many 
variables such that analytical solutions are difficult to obtain and numerical solutions are 
prone to errors. Further, many PTM systems admit multistationarity (the existence of 
more than one steady state under particular biological conditions) which is a mechanism 
for cellular decision making [18j . It is therefore of interest to determine the parameters for 
which mono- and multistationarity occur. Several non-necessary conditions for a unique 
positive steady state are known [H 131 El) but when these fail, multistationarity is difficult to 
determine and often decided based on a random parameter search. Procedures to eliminate 
variables (hence, equations) is therefore fundamental to the theoretical understanding of 
these systems as well as for numerical analysis. 

Our work is inspired by previous work by Thomson and Gunawardena (TG) [17J which we 
extend to embrace a range of important PTM systems such as signaling cascades (including 
the MAPK cascade) and two-component systems with phosphorelays and phosphotranfer 
|14| . as well as systems with self- interactions. We develop the idea of a cut Sa, a subset of 
the substrates S with certain properties that allow us to express the steady state equations 
as rational functions in the "core" variables S\Sa, providing an algebraic parameterization 
of the steady states in terms of the core variables. If the core variables take positive values 
at steady state, then we show that all other concentrations are either zero or positive as 
well. 

Further, we show that cuts relate to conservation laws (conserved quantities that imply 
that the dynamics takes place in an afhne invariant subspace of M") that arise as connected 
components in the species graph |lj. Conservation laws are often used as a first step to 
reduce the dimensionality of the system. In our approach, conservation laws come into 
play after elimination of variables from the steady state equations. In this way, we allow 
for a larger reduction in the number of core variables. 

Our appoach makes use of algebraic tools as well as some basic graph properties; for 
example Tutte's Matrix- Tree theorem |19l I17j . One benefit is that parameters are treated 
as symbolic constants and do not need to be fixed or assumed known in advance. This 
is particularly relevant in biology, where we often are faced with systems that depend on 
experimental parameters (kinetic rates), which are difficult to determine. 

2. Post-translational modification systems 

2.1. PTM system. A post-translation modification (PTM) system consists of two non- 
empty sets of species, S = {Si, . . . , Sn} (the substrates) and y = {Yi, . . . , Yp} (the inter- 
mediate complexes) with 5 n 3^ = 0, and a set of reactions Ret = RaU RtU RcU R^ with 
associated positive reaction rate constants: 

Ra = {Si + Sj ^ Yk\{i, j, k) G la} Re = {Yi ^ Yj I (i, j) elc,i^ j} 

i.k 

Rb = {Yk ^S^ + Sj I {i, j, k) G h} Rd = {Si ^ Sj \{i,j)eld,i^ j} 

ioT laJb ^ {l,...,iV}' X {!,..., P}, /, C {!,..., P}2 and/d C {!,..., iV}2. To fix the 
notation, we assume that any {i,j,k) £ Iq U /^ satisfies i < j, so that self-interactions a 
priori are allowed. If the rate constants are not required, we put an arrow to indicate a 
reaction and omit the rates. Further: 
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(i) All chemical species are involved in at least one reaction. 

(ii) For every intermediate complex Y^ there exist i < j, indices ki, . . . ,kr and a chain of 
reactions 1^ — )• 1^^ —)■•••—)• Yi-v -^ Si + Sj . 

Assumption (ii) ensures that Y^ ultimately dissociates into two substrates. Also, we allow 
that there are more than one Y^ such that Si + Sj — )• 1^ or Y^ ^ Si + Sj for given Si , Sj . 
For convenience, we put Cij = 0, dij = if (i,j) ^ Ic or I^ respectively, and similarly 
a^j = and 6^ • = if {i,j, k) ^ la or /(,, respectively. For i < j and k, we define a^^ = af • 
and 6^,- = b^ ,. For later use, we define 

5(3 = {Si G 5| (i, i, k) £ laU lb for some A:} 

to be the set of self-interacting substrates. 

This setting fits post-translational modification of proteins catalyzed by enzymes as well 
as the transfer of modifier groups: 

E + S^^Y ^-E + S* P* + S^^Y^^P + S* 

where S* ,P* are modified proteins (substrates), S, P their corresponding unmodified forms, 
E an enzyme (substrate) and Y an intermediate complex. That is, 5 = {S,S* ,P,P* ,E} 
and y = {Y}. In the first case, the attachment of the modifier group is catalyzed by 
the enzyme E, whereas in the second case, a modifier group is transferred from P* to S. 
Modification of a substrate or an intermediate complex without the involvement of other 
species is modeled by 5 ^ 5* and Y ^ Y* , respectively. 

As an example consider the PTM system with S = {Si, S2, Ss, ^4, 55}, y = {Yi, Y2, Is} 
and reactions 

I tin o ("i o /> 

(2.1) Si^S2 S2 + S3^=^Yi:^=±Y2^^Si + S4 

S4 + S5 ^^ Y3 —^ S3 + S5 

One interpretation is that Si is modified to ^2. The modifier group is then transferred from 
5*2 to 5*3 with the formation of two intermediate complexes Yi,Y2, causing the modification 
of 5*3 to 5*4 and the demodification of 5*2 to Si. Finally, 54 is demodified via a Michaelis- 
Menten mechanism catalyzed by an enzyme 55. 

Nomenclature. We introduce a few concepts that will be used in the following, some 
of which are taken from Chemical Reaction Network Theory (CRNT) [6j[8j. Consider the 
set of complexes of the reaction system: 

C = yu {Si, Sj\ {i,j) e Id} U {Si + Sj\ {i,j, k)eIaU h for some k). 

Then: 

• A £ C reacts to B £ C if there exists a reaction A ^ B. 

• A £ C ultimately reacts to B £ C if there exists a sequence of reactions j4 —)• ^1 —)•••• —>■ 
Ar ^ B with Am £ C. If Am £ y '^y for all m, then A ultimately reacts to B via y. 

• Si and 5,- interact if for some Yk either Si + Sj reacts to Yfc or vice versa. 
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• Si,Sj are 1-linked if dij or dj^i ^ 0. Yf:,Yy are 1-linked if Ck^v or c„^fc / 0. Si and Y^ are 
1-linked if for some j, Si + Sj reacts to Y^ or vice versa (j = i is allowed) . 

Assumption (ii) of a PTM system ensures that all intermediate complexes ultimately 
react to some Si + Sj via 3^. 

2.2. Mass-action kinetics. The set of reactions together with their associated rate con- 
stants give rise to a polynomial system of ordinary differential equations taken with mass- 
action kinetics: 

N j P 

Yk = Yl ^{aljS.S, - 6f,,n) + ^ic,,kY, - Ck,,Yk), k = l,...,P, 

j=l i=l ti=l 

N P N 

Si = Y,Yl ^^'j- (-4jSiSj + ^tjYk) + Y.(^j^^^j - '^i'A)^ z = 1, . . . , TV, 

3=1 fc=l j=l 

where e^j- = 1 if i ^ j and 2 if i = j and where x denotes dx/dt for x = x{t). Here we abuse 
notation and let Si^Y^ denote the concentrations of the species Si^Y^ as well. The steady 
states are the solutions to the polynomial system obtained by setting the derivatives to 
zero, i.e. Y^ = and Si = 0: 

N j P 

(2.2) = Y, Y.^(4,J^iSj - hl^Yk) + Y,{c,,kY. - Cfc,,yfc), A; = 1, . . . , P, 

j = l i=l v=l 

N P N 

(2.3) 0=Y,Y1 ^iA-4jSiSj + tj^k) + Y.{dJ,^S, - di^,Si), i = l,...,N. 

j=l k=l 3 = 1 

This system is quadratic in the variables Y^^Si, but the only quadratic terms are of the 
form SiSj. It is linear in Yk- 

It is convenient to treat the reaction rate constants as parameters with unspecified 
(positive) values and view a^ ,,h\^,Cu,v,dw,t as symbols. For that, let 

Con = {a\ Mi, j,k) G /J U {6^.,|(i,i, /c) G /&} U {cfc,^|(A;,i;) G /J U {4,^|(A;,^) G Id} 



be the set of the non-zero parameters (symbols). Then, the system (2.2)-([2^ is quadratic 
in 5 U 3^ with coefficients in the field M(Con). Further, if all Si are considered part of 
the coefficient field, then the system is linear with coefficients in M(ConU5) and variables 
Yi,...,Yp. 

Only non-negative solutions of the steady state equations are biologically meaningful. To 
study positivity of solutions, we introduce the concept of S-positivity. Let X = {xi, . . . , Xr} 
be a finite set. A non-zero polynomial in M[X] with non-negative coefficients is called S- 
positive. Similarly, a rational function / is S-positive if it is a quotient of two S-positive 
polynomials. If xi, . . . , x,- are substituted by positive real numbers in /, we obtain a positive 
real number. In general, a rational function / = p/q in zi, . . . ,Zs and coefficients in M(X) 
is S-positive if the coefficients of p and q are S-positive rational functions in xi , . . . , x^ . If / 
is a rational function in xi, . . . , x,. and Xj = g{xi, . . . , Xj, . . . , Xr) with g a rational function, 
then substituting g into / gives / as a rational function in xi, . . . , Xj, . . . , x^. 



The differential equations of Example (2.1) are: 
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(2.4) Yi = al^S2S:i - (&^,3 + 01,2)1^1 + C2,iY2 Si = -di,2Si + bl^Y2 

Y2 = Cl,2>l - (fc?,4 + C2,l)l2 ^2 = di^2Si - al^^S2S3 + ^a.s'^^l 

n = al.S^S^ - (61,5 + ^3,5)n ^3 = -al3S2S3 + bl^Yi + bl.Y^ 

S^ = -al^SiS^ + (6l,5 + 61,5)^3 ^4 = -al,S:,S5 + 6?,4y2 + ftls^s- 

To compute the steady states, we can use I3 = to eliminate Ya as a function of the 
substrates. Also Yi,Y2 can be eliminated by solving the linear system Yi = Y2 = 0. This 



is a general feature of PTM systems and is covered in Section [3^1 

Further, observe that 55 + ^3 = 0, which implies that the sum ^5 + I3 is independent of 
time and thus conserved. In fact, it implies that one of the equations Y^ = and 8^ = 
is redundant. Removing one of them leaves a polynomial system with 7 equations in 8 
variables, and thus the solutions to the steady state equations form an algebraic variety 
of dimension at least one. This redundancy can be compensated for by fixing the value 
S5 + Y3 = S and adding this relation to the steady state equations. 

In the next section we discuss the existence of the so-called conservation laws and provide 
a graphical procedure to determine (some of) them. In most cases the procedure provides 
a set of independent conservation laws, but, as will be discussed below, this might not 
always be the case. 

2.3. Conservation laws. We consider systems where inflow of species is not allowed and 
species are not degraded or able to diffuse out. Such systems are "entrapped" in contrast 
to open systems (so-called "continuous flow stirred tank reactors") [3]. PTM systems are 
entrapped and have conservation laws that reflect that the total amount of species remains 
constant either in free form Si or in bounded form Yj. These laws follow from the system 
of differential equations and appear as linear combinations of species (e.g. ^5 + Is = S" in 
the example above). 

The existence of conservation laws implies that the dynamics of the system takes place 
in a proper invariant subspace of M^"*"^. We identify M.^'^^ with the real vector space 
generated by 5 U 3^ so that M^+^ = {Si, . . . , Sn, Yi, . . . , Yp). The species Si and Yk are 
unit vectors with a one in the i-th and (A^ + fc)-th entry, respectively, and all other entries 
being zero. A vector v = {Xi, . . . , Xn, fJ-i, ■ ■ ■ , fj.p) is identifled with the linear combination 
of species Y^^ XiSi + Y^k t^kYk- 

Consider the stoichiometric subspace of M "*" [3] of a PTM system: 

r = {Si + Sj - Yk\ ii,j, k)eIaU h) + {Yk - y.l {k, v) e h) + {Si - 5,| (ij) g Id). 

If (Ai, . . . , Xn, /ii, • • • , /Up) G P-*-, then ^- XiSi + ^^^ /ifelfc = 0. The converse might not 
be true [8]. It follows that any basis {uj^, . . . ,(^'^} of P"*" provides a set of independent 
conserved quantities ^j=i A-S^ -|- X]fc=i fJ-i^k if "^^ = (Ai, . . . , A^y, M, • • • , Aip)- Therefore, 
if total amounts Si, . . . , Sd G K+ are provided, we require the steady state solutions to 
satisfy: 

N P 



(2.5) Si = J2^'iS, + Y,fiiYk 1 = 1,. ..,d. 



k=l 
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Total amounts are fixed by the initial concentrations of the species. We say that equations 



(2.5) are independent if the system has maximal rank, or equivalently, if the corresponding 
vectors of F-*- are independent. 

We introduce the concepts of a cut and a non-interacting graph and show that they 
provide means to obtain conservation laws. 

Definition 2.6. Let a non-empty set 5q, C 5 be given and let the associated set 3^^ C y 
be the smallest set such that Yk E 3^^ if Ifc is 1-linked to some 5^ G Sa or to Ym G 3^a- 

(i) Sa is dosed if Sj belongs to Sa whenever Si £ Sa is 1-linked to Sj, and if Si and 

Sj interact and are 1-linked to Y^ G 3^a, then Si or Sj are in Sa- 
(ii) Sa is a cut if (a) Si, Sj G Sa do not interact for any i,j, and (b) Sa is closed, 
(iii) A cut Sa is minimal if it has no proper closed subsets. 

Condition (ii) implies that a self-interacting substrate S £ Sq cannot belong to any cut, 
that is, 5a n 5(3 = for any cut Sa- Note that a closed subset S' of a cut is also a cut. 
The union of two disjoint cuts Sa,S'^ is a cut if 3^0 n 3^^ = 0. 

In the PTM system with reactions Si + S4 ^s^ Y2 ^^ S2 + 84, and Yi ^5^ ^2 + ^3 , the 
set {<S'i,52} is a cut, while {S'ljS's} is not. There are no proper closed subsets of {5i,S'2} 
and thus the cut is minimal. 

Definition 2.7. Let a non-empty set 5q C 5 be given and let 3^0 C 3^ be as in Defi- 
nition 2.6 Further, let Gs^y^ be the graph with node set Sa U 3^^ and edges between 



1-linked nodes. The graph is non-interacting if it is connected and Sa is a cut. 

If Sa = S, then ya = y. All graphs Gs^y^ are naturally subgraphs of Gsy- Without 
proof we state the following: 

Lemma 2.8. Let Sa be a cut and G' be a connected subgraph of Gs^y^ with node set 
S' U y , S' C Sa and y' C ya- The following are equivalent: 

(i) S' is closed with associated set y' . 
(ii) G' is a connected component ofGs^y^- 
(iii) G' is non-interacting and contains only species in Sa U 3^q. 

// either is the case, then S' is a minimal cut and G' = Gs'y - 

Thus, the non-interacting graphs containing substrates only in a cut Sa are exactly the 
connected components of Gs^y^- All non- interacting graphs contain some node from S 
(condition (ii) of a PTM system). However, such a graph might not exist. Consider for 
example the system with reactions ^i ^^ 53, 82^ S'^, Si + S2^Yi. The graph G^^y 



IS 



Si 
Yi^ ^S, 
82^ 

Condition (b) of Definition |2.6[ ii) implies that any non- interacting graph must contain all 
four species, which contradicts condition (a) of the same definition. 



7 Variable elimination in post-translational modification systems 

Lemma 2.9. Let Hi, . . . , Hn be the non-interacting graphs of a PTM system, Ci the node 
set of Hi, Si = S nCi and 3^; = y CiCi. Then, a)/ = for 

ioi = Y,s+Y,y 1 = 1,. ..,n. 

SeSi YeVi 
That is, Hi corresponds to a conservation law and uji is fixed by the initial amounts. 

Proof. Substrates in C; interact only with substrates in S \ Si and thus, by definition of 
yi, if a^ • 7^ or 6f ,• 7^ for i y^ j then: (a) if Si (resp. Sj) is in Si, then Sj (resp. Si) 
belongs to S \ Si, and Y^ £ yi; (b) if Yk £ yi, then either Si or Sj , but not both, belongs 
to Si. If c^_fc / or Ck,v 7^ 0, then Yk, Yy belong to the same non-interacting graph (if any); 
if dij / or dj^i ^ 0, then Si, Sj belong to the same non-interacting graph (if any). Since 
5/ n 5(3 = for Yk £ yi and Si G Si we have: 

yk= Yl YI i4j^i^j-''idYk)+ YI KkY.-ck,vYk) 

Si= Y E {-4jSiSj + b'^^^Yk)+ Y (dj^iSj-d^jSi). 
k\Yk&yi j\Sjes\Si jls.eSi 

It follows that EklY^ey^ Ev\Y^GyMv,kYy-Ck,vYk) = and Ei|5.e5, J2j\s,eSiidj,iSj-dijSi) = 
0. Similarly, the remaining terms in ui cancel. Thus, oJi = 0. D 

Thus, each non- interacting graph gives rise to a conserved amount. If each non-interacting 
graph contains a species that only belongs to that specific graph, then the w^'s are indepen- 
dent. In particular, conservation laws derived from the connected components of Gs^y^ 
for some cut Sa are independent. In general, the set of conservation laws found from 
Lemma |2.9| can be reduced to a set of independent conservation laws. 



In Example (2.1), the graph Gsy is 

H3 




The non-interacting graphs Hi, H2, H3 are colored. If total amounts Si, S2, Ss are provided 
then the steady state solutions must satisfy: 5i = 5*5 + 13, S2 = Si + S2 + Yi + Y2, 
and 53 = S3 + S/i + Yi + Y2 + Y^. These conserved total amounts are easily verified by 



differentiation using (2.4). 

Consider a two-layer cascade of modification cycles that share the demodification en- 
zyme F in each layer. The reaction system consists of 5 = {E, F, Si, S2, S3, S4}, y = 
{Yi,Y2,Y3,Y4} and the reactions 

(2.10) E + Si^^Yi^E + S2 F + S2:^Y3^F + Si 

S2 + S3^±Y2^S2 + Si F + Si^±Y4^^F + S3 

The subsets Sa = {E, S3, S^}, {E, F}, {Si, S2} are examples of maximal cuts (they cannot 
be extended to larger cuts by including more substrates). The graph Gs,y is 
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These graphs are obtained as connected components of the graph Gs^y^ for the cuts Sa 
above. As in the previous example, the different non-interacting graphs yield independent 
conservation laws and thus if total amounts are provided, we obtain the following equations: 
Si = E + YuS2 = S3 + Si + Y2 + Yi,S3 = F + Y3 + Yi and S^ = Si + S2 + Y1 + Y2 + Y3. 
This procedure provides an easy construction of conservation laws. In the two examples 
above, the conservation laws obtained from the graph are independent and, additionally, 
determine all conservation laws arising from T-^ (dimr-*- = 3 and 4, respectively). However, 
this is not always the case. Consider for instance the reaction system 



(2.11) 

The graph Gsy is 



Si + S2:^Yi^^Y2^±S3 + S4 

^ yi - ^2 ^ 

'S'2 S3 

There are 4 non-interacting graphs that give the conserved total amounts Si = Si + S3 + 
Yi + Y2, S2 = Si + Si + Yi + Y2, S3 = S2 + S3 + Yi+ Y2, and Si = S2 + Si + Yi + Y2. The 
rank of the space generated by the corresponding 4 vectors in M® is 3, implying that one 
of the relations is redundant. In this case the procedure still gives all conservation laws, 
because the dimension of F is 3. 

Consider the following reaction system: 



(2.12) S1 + S3 

The graph Gs,y is 



Yi^±S2 + S4 



Si + 84^^ Y2 



52 + 53 ^ Y3 



'S'l 'S'2 

^2 \ / ^1 ^ / ^3 

There are 2 non-interacting graphs that give the conserved total amounts 5i = S'l + 
S2 + Y1 + Y2 + Y3, and S2 = S3 + Si + Yi + Y2 + Y3. However, dimT^ = 3 and the 
procedure fails to provide three independent conservation laws. A third conservation law 
is 53 = S'i-|-S'4-|-yi-f2y2j and the coefficient 2 of 1^2 cannot be obtained from non-interacting 
graphs. 

2.4. Cuts of S and conservation la'ws. We provide a criterion to guarantee that there 
are dimr-*- independent conservation laws derived from non-interacting graphs. The crite- 
rion will be used in Section [H 

In the following we make use of Lemma 2.8 without further reference. Let Sa be a cut 
with associated set 3^^. Define S^ = S\Sa and y^ = y\ya, and let Na, Pa (resp. Nl^, P^) 



be the cardinality of Sa, 3^a (resp. 5^, 3^^). Extend the set of conservation laws derived 
from the connected components of Gs^ y^ to a maximal set of n independent conservation 
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laws derived from other non-interacting graphs (thus containing species in S^ U y^)- Let 
n'^ = n — Ua, where tIq, is the number of connected components of Gsa,ya- 

Lemma 2.13. Let Sa be a cut and keep the notation introduced above. Then, we have that 
dim {{S^ U y^) n r) < A''^ + -P^ — "•« and dimF-*- = n if and only if 

dim ((5^ u 3^^) nr) = N^ + p^- <. 

Proof. Without loss of generality we can assume that ya = {Yi, . . . ,Yp^} and Sa = 
{Si,...,SnJ. Identify R^+P with M^« x M^« x M^S x M^^ and let 

Ta = {A- B\ for each edge A — B in Gs^^^). 

The space T-^ is generated by the vectors which are sums of species in each connected 
component oiGs^y^ and hence dimF^ = na- We have dimr-*- > n = na+n% and we want 
to determine when equality holds. Equivalently, we want to see when dim T = N + P — n. 
If this is not the case, then dim T < N + P-n. Note that N = Na + N^ and P = Pa + Pa- 
Note that dim Ta = Na + Pa- ria- Let vr : M^+^ -^ ^No,+Po, denote the projection onto 
the first Na + Pa coordinates and 7r„: F — )• F^ its restriction to F (tTq, a surjective map). 
Then, dimF = dim F^ + dim ker tTq and so dimkervTo, < N^ + P^ — n^. Further, dimF-*- = n 
if and only if dim ker tTq = A''^ + P^ — n^ . Finally, note that (5^ U J'^) fl F = ker tTq, . Indeed, 
let z: F <^ M^+^ and ia'- ^a ^^ M^"+-^" denote the natural inclusions. We have that 
ia° T^a = vr o i. The kernel of vr is clearly M 'S+^a = (^^^ u y^'j from where it follows that 
the kernel of iTa is (5^ U 3^^^) n F. 

Therefore, dim ((5^ U 3^^) n F) = dim ker vTq, = A'^^ + P^ — n^ if and only if dimF"*- = n 
and the lemma is proved. D 

As each non- interacting graph corresponds to a minimal cut, the lemma above provides 
a condition for when all conservation laws are recovered from cuts. 

Remark. An easy way to construct elements of {S^ U 3^^) n F is by considering: 

(i) Vectors Si — Sj for any pair Si , Sj G S^ for which there exists a chain of reactions 

Sm + Si — Ai — • • • — Ar — Sm + Sj for some Sm, Au £ C and — is -^ or — )•. 
(ii) Vectors Si + Sj — Yfc, Si — Sj or Y^ — Y^ corresponding to reactions with Si, Sj G S^ 
andyfc,^; e3^a- 
If we can construct A^^ + P^ — n^ independent elements of (5^ U 3^^) n F of the previous 
type, then the previous lemma holds. 



In Example (2.1) consider the cut Sa = {Si, 5*2,55} with S^ = {53,54} and the given 
conservation laws (n = 3). We have A"^ = 2 and n^ = 1. Further, ya = y so that P^ = 0. 
The element 53 - 54 = (53 + S5 - Y3) - {S4 + S5 - Y3) belongs to (5^) n F. In addition, 
^a + -fa ~ ^a = 1 ^''^^ ^hus diin{{Sa) Pi F) = A''^ + P^ — n'^, implying that all conservation 
laws are found from non- interacting graphs. 



In Example (2.10), consider the cut Sa = {E, S3, S4} with Sa = {F, Si, S2} and A"^ = 3. 
In this case, 3^^ = {^1)^2,^4}; 3^q = {^3} and so P^ = 1. Two of the four conservation 
laws involve elements in Sa U ya only and hence Ua = na = 2. Further, A"^ + P^ — n^ = 2. 
The two independent vectors P + 52 — V3 and P + 5i — Ys belong to (P, Si, S2, Y3) n F. 
Thus, the graphical procedure provides all conservation laws. 
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In Example (2.11), consider the cut Sa = {81,83} with S^ = {52,54} so that N^^ = 2, 
P^ = 0. There is only one conservation law in 5q, U 3^^, 53 + 51 + 11 + 12, and since n = 3, 
then n^ = 2. It follows that A'^ + P^ — n^ = 0, and we are guaranteed that the dimension 
of (52, 54) n r is zero. 



In Example (2.12), consider the cut Sa = {5i,52} with 5^ = {53,54} and N^^ = 2, 
P^ = 0. We have n^ = 1 and N^ + Pa~''^a — 1- However, (53, 54) n F has dimension zero 
and thus not all conservation laws arise from non- inter acting graphs. 

3. Variable elimination 
In this section we show that the intermediate complexes can always be eliminated and 



expressed as polynomials in the substrates with coefficients in ]R(Con) (Section 3.1). After 



choosing a cut Sa , the substrates in Sa can be expressed in terms of those in S^ = S \Sa 



(Section 3.3) 



3.1. Elimination of intermediate complexes. Consider the system 1^ = in ( |2.2[ ) as 
a linear system of P polynomial equations with coefficients in M[ConU5] and P variables 
Yi, . . . , Yp. If the system has maximal rank, then there is a unique solution in M(Con U S). 
Specifically, we have a linear system AY = z where Y = (Yi, . . . ,YpY and A = {Xk,v} 
is a P X P matrix with coefficients in M[Con], 

_ j-Cv,k if k^v 

'""lEiLiELife^ + E^Lic,, ifk = v. 

The independent term z = (zi, . . . , zp)* is in M[ConU5]: Zk = Yli<i '^ij^i^j- 

Assume that A has maximal rank P in M(Con). Then, using Cramer's rule to solve 
linear systems of equations, we obtain that 1^ = pk/p with p = det{A) 7^ and pk 
the determinant of A with the k-th. column substituted by z. Since the determinant is 
a homogeneous polynomial in the entries of the matrix, it follows that p £ M[Con] and 
Pk G M[ConU5]. Therefore, 

Yk = ^2__^^j jOjOj 

with fi^j G M(Con) and thus Yk is a polynomial in M(Con)[5]. If both p, p^ are 5-positive 
elements of M[Con] and M[ConU5], respectively, then for positive rate constants and non- 
negative values of 5j, the steady state value of Yk is non- negative as well. S-positivity of 
p, (T is proven in the next section using the Matrix-Tree theorem [19]. Some basic concepts 
from graph theory are required. 

Graphs and the Matrix- Tree theorem. Given a directed graph G, a spanning tree 
T is a directed subgraph with the same node set as G and such that the corresponding 
undirected graph is connected and acyclic. There is a unique undirected path between any 
two nodes in a spanning tree [5] . A spanning tree r is said to be rooted at a node v if the 
unique path between any node w and v is directed from w to v. It follows that v is the 
only node with no out-edges, that is, there is no edge of the form v ^ w in t. In addition, 
there cannot be a node with two out-edges in r. The graph G is strongly connected if for 
any pair of nodes v, w there is a directed path from v to w. Any directed path from v to 
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w in a strongly connected graph can be extended to a spanning tree rooted at w. Some 
general references for graph theory are [5J and |10j . 

If G is labeled, then r inherits a labeling from G and we define 

^(^) = n "• 

a 
X — tydT 

Assume that G has no self- loops. Order the node set {vi, . . . , Vn} of G and denote by atj 
the label of the edge Vi — ;■ Vj. We set aij = if there is no edge from Vi to Vj (thus 
0'i,i = 0). Let JC.{G) = {ctij} be the Laplacian of G, that is the matrix with 



a 



*j 






such that the column sums are zero. For each node Vj, let Q{vj) be the set of spanning 
trees of G rooted at Vj. Then, the Matrix- Tree theorem states that the maximal minor 
C{G)uj\ (the determinant of the minor obtained by removing the i-th row and the j'-th 
column of C(G)) is: 

£(G)(,,.) = (-l)""i+^+^- Y. ^(^)- 

Note that for notational simplicity we have defined the Laplacian as the transpose of how 
it is usually defined and the Matrix- Tree theorem has been adapted consequently. 

In our case, the matrix A is not a Laplacian, since the column sums Ylij=i Yll=i ^ij ^^^ 
not zero. However, A can be extended such that its determinant is a maximal minor of a 
Laplacian. 

3.2. Decomposition of the system. Let Gy be the directed graph with node set y and 
a directed edge Y^ — )• Yy if {k,v) £ Ic- The node sets of the connected components of Gy 
determine a partition of 3^: 3^ = 3^i U • • • U 3^^- Let Pi be the cardinality of 3^/ and rename 
the intermediate complexes such that yi = {Yp-^^ |_Pj_^+i, . . . , Yp^-| \-Pi}- 

If Yk G yi for some /, then Ck,v = Cy^k = for any v such that Y^ ^ yi. It follows, that 
A is a block diagonal matrix diag(^i, . . . ,As) with Ai being a Pi x Pi matrix. Solving 
AY = z is thus equivalent to solving s "smaller" systems with matrices Ai. Further, A has 
maximal rank P if and only if Ai has maximal rank Pi for all I. 

Consider the connected component Gy^ corresponding to 3^;. We construct an extended 
labeled directed graph Gy^ with node set yi U {*}. For convenience we order the nodes 
such that Yp^-f ^Pj_^^k is the k-th. node and * the {Pi + l)-th node. Let b^ = Yli<j ^i,j 

and a = z_ji<j '^i j- The graph Gy^ has the following labeled directed edges: Y^ — ^ Y^ if 
(A;, v) G Ic, Yfc -^ * if ft V 0, and * A Y^ if a^ / 0. 



In Example (2.1), the graph Gy has two connected components 3^1 = Yi < ^ Y2 and 



y2 = Y3. The graphs Gy-^ and Gy^ are 
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Cl,2 




Y-, 



''3,5+''4.5 



554S5 



Let £ = {ok^v} be the Laplacian of Gy, . If k,v < Pi, then a^,?; = —^k,v The entries of 
the last row are ap^+i^fc = b^ for k < Pi and the entries of the last column are a^^Pi+i = a^ 
(= Zk) for k < Pi- We conclude that the (P; + 1,P; + 1) principal minor of C is exactly 
—Ai and thus, by the Matrix- Tree theorem, we have 

det(AO = (-l)^'£(P,+i,P,+i) = Yl <^)- 

ree(*) 

Assumption (ii) of a PTM system ensures that each Y^ ultimately reacts to some Si + 5^ 
via y, and hence there exists at least one spanning tree rooted at *. Thus, det(^j) 7^ 
and (lei{Ai) is an S-positive element of M[Con]. 

By the definition of pk and the Matrix- Tree theorem, 

pu = i-l)'+'C^p,+i,k) = E ^(^)' 

ree(yp,_j+fe) 

and hence pk is either zero or an S-positive element of M [Con U 5]. 

If there exists at least one spanning tree rooted at Vk = yp,_j+fc, then pfc / 0. A necessary 
condition for this to happen is the existence of at least one in-edge to Vk- Otherwise the 
concentration at steady state of v^ is zero, which is expected if Vk is only consumed and 
never produced. Similarly, if there is no reaction of the form Si + Sj — )■ Yp^_-^^rn for any m 
(that is, a directed edge * — >• Vm), then pk = for all k. 

The term p^ is a homogeneous polynomial of degree 2 in 5 with coefficients in ]R[Con], 
because any spanning tree rooted at a node Vk has exactly one edge of the form * — >• f m 
for some m. Further, a monomial SiSj appears in pf^ only if Si + Sj ultimately reacts to Vk 
via yi. If Gy, is strongly connected, then this condition is both sufficient and necessary. 
Indeed, if 5^ + Sj ultimately reacts to Vk via yi, then there is a spanning tree rooted at Vk 
containing this path. 

The next proposition summarizes the discussion above: 

Proposition 3.1. Consider a PTM system with intermediate complexes y and substrates 
S. Then, Y^ = for all k, if and only if 

(3.2) Yu = Y,pIjS,S, 



with pf , G M(Con) being either zero or S-positive. Further: 



ij 



If Si + Sj does not ultimately react to Y^ via y, then p^- = 0. 



(ii) If Gyi is strongly connected andY^ G 3^;, then p^ ■ / if and only if Si + Sj ultimately 
reacts to Y^ via yi . 



(Hi) Gyi is strongly connected if and only if in (3.2), Y^ is a non-zero polynomial in 
M(Con)[5] for all Yk £yi. 
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Remark. The condition that Gy^ is strongly connected is biochemically reasonable: 
The intermediate complexes are not the initial or final products of the system and should 
eventually be broken up into parts. 

In Example (2.1), the graph Gy-^ has three spanning trees rooted at * so that det(j4i) = 



bi 4Ci,2+^2 3^2,1+^1 4^2 3- There is one spanning tree rooted at Y2, giving p2 = ci,202 3'5'25'3, 
and two spanning trees rooted at Yi, giving pi = (6f 4 + £2,1)02 3'S'25'3. The graph Gy^ has 
one spanning tree rooted at * so that det(742) = 635 + 545, and one spanning tree rooted 
at Y3, giving p^ = al^^SiS^. Thus: 

^1 = /^2,3'5'25'3, Y2 = ^2,3^283, Y^ = p^^^SiS^ 

Wixn ^2,3 - det(Ai) ' ^2,3 - det(Ai)' ^"^^ /^4,5 - det(A2)- 

Lemma 3.3. Let Gy = ^iGyj^ . The graphs Gy^, I = 1, . . . ,s, are strongly connected if and 
only if the graph Gy is. 

Proof. Assume that the graphs Gy^ are strongly connected. Then, for any v G Gy^ and 
uj £ Gy^ , there are directed paths i; — )• * in Gy, and * — ;■ u in Gy. , which by composition 
give a directed path between v and w. 

For the reverse implication, let v,uj be two elements of yi. Since Gy is strongly con- 
nected, there exists a directed path a: v ^ w in Gy. We can assume that v,w ^ *. A 
path connecting an intermediate complex in 3^/ to one in yj for j 7^ / must pass through 
*. If a path a goes through v £ yj, for j ^ I, then it must go through *, first in and 
then out, potentially many times until it goes back to yi and to w. Therefore, a has the 

form V —h * — > * — )• w with ai and 02 being paths in Gy^. It follows that the path 
V — > * — > w is a. directed path from t> to w in Gy, . D 



3.3. Elimination of substrates. Equation (3.2) shows that at steady state the interme- 



diate complexes are given as zero or S-positive rational functions in the substrates and the 



rate constants. Insertion of (3.2) into the (time dependent) differential equations for the 
substrates is the procedure known as the quasi-steady state assumption. The rationale is 
that intermediate complexes tend to reach steady state much faster than substrates and 
thus some variables in the dynamical system can be eliminated. We have shown here that 
PTM systems "mathematically" enable this simplification although justification is required 
in concrete examples. 



We now use the steady state equation (2.3) to further eliminate some of the substrates 



in terms of others. Recall equation (2.3), that is 5^ = 0, 

N P N 

(3.4) = E E 'iA-4JS^Sj + b'^jYk) + Y,{dj,^Sj - di,,Si) 

j=i k=i j=i 

for i = 1, . . . ,N. After substitution of the values for Y^, we have 

TV P N P N 

(3.5) o = Y,Y.Yl ^^Aui^itSjSt - E E ^^A,s^sJ + Y.(dj,Sj - d,,,s,), 

u=l fc=l j<t j=l fc=l j=l 



14 



Variable elimination in post-translational modification systems 



These equations are quadratic in 5. To proceed with hnear ehniination it is necessary to 
decide which variables are to be eliminated and which will be taken as part of the coefficient 
field. Since a monomial SiSj appears only if Si and Sj interact, we can proceed as long as 
S can be partitioned in an appropriate way. 



Lemma 3.6. Assume that there is a cutSa with associated sety^ (Definition 2.6). Then, 
(3.5) for the substrates in Sa is a homogeneous linear system of equations in the substrates 



Sa and with coefficients in M(ConU5^). 

Proof. For the three sums in ( |3.5[ ) we make the following observations: If Si E Sa and 
dj.i / 0, then also Sj G Sa. If Si G Sa and of,- / 0, then Sj Sa, otherwise Si and Sj 



-^j,i 



would interact. Finally, if n^ 



■j,t 



/ 0, then according to Proposition 



3.1 



Sj + St ultimately 

reacts to Yfc via 3^q,. Hence, since Sa is a cut, one of Sj and St (but not both) belongs to 
Sa. Thus (3.5) for the substrates in Sa is a homogeneous linear system of equations in the 
species in 5^. D 

Assume that there exists a cut Sa and that Sa = {Si, . . . ,S]\f^}. It follows that for 



Si G Sa the equations in (3.5) form an Na x Na homogeneous linear system of equations 
with variables Sa and coefficients in M(ConU5^). Further, eij = 1 if Si £ Sa- 

Oj for i = j, where 



Let B be the matrix with entries bij for i ^ j and b, 

Ni N 



E^M+ E 



p 



Sj , b 



*j 



d 



'3,1 



j=Ni+lk=l 



N P 
t=Ni+l fc=l 



TV 

E 

u=Ni+l 






so that (3.5) becomes 



(3.7) 



E KjSj^{h, 



ai)Si. 



Consider Example (2.1) and the cut Sa = {Si, S2, S^}. Then the equations (3.5) are 
= — (ii,2'S'i + bi ^iJ-2,s>^2S3 and = ^1,25*1 — a2 3'5'2'S'3 + ^2 3/^2 3'S'2'S'3, corresponding to 
Si = and 5*2 = 0, respectively. The equation 5*5 = is trivial because of the conservation 



law Y3 + S5 = 0. Further, we have ai = di 2, 02 



''2.3 



5*3, bi,2 — &l,4/^2,3'S'3 5 ^2,1 — C?l,2, 



^2,2 = ^2 3/^2 3^3^ while the rest of the coefficients are zero. 



Lemma 2.8 ensures that there is a conservation law for each connected component of 



Gsa,ya- Let Cf , . . . , C^^ be the node sets of the connected components and define Sa^i = 
Sa n Q° and ya,i = yaD Cf so that Y^SieS^, i ^i + Y^Y^^y^ ^ Yfc = for / = 1, . . . , n^, are 
conservation laws. Imposing only that the intermediate complexes are at steady state, that 
is y^; = for all A;, we obtain 



(3.8) 






5. = 0, 



1 = 1, 



,nc 



It follows that the column sums of the matrix B restricted to the rows corresponding to 
the substrates in Sa,i are all zero. Consequently, the matrix B has rank at most Na — n^. 
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L et G y^ be Gy restricted to the nodes ya- It follows from the definition of 3^^ (Defini- 
tion 
Lemma 



2.6) that Gy^ is a union of connected components of Gy. Define Gy^ similarly (cf. 



3.3). Let Na^i be the cardinality of Sa^i- 



Lemma 3.9. After reordering the substrates in Sa, B is a block diagonal matrix, namely 
diag(i?i, . . . ,Bna), where Bi is an N^^i x N^^i matrix. Further, if bij ^ then there is a 
reaction Sj — )■ Si or there exist Su, St £ S^, so that Sj + St ultimately reacts to Si + Su via 
ya. If in addition Gy^ is strongly connected, then the reverse is true. 



Proof. It follows from Lemma 3.3 and Proposition 3.1 h) that if /i^^ 7^ then Sj + St 
ultimately reacts to Y^. By definition, b^ ^ ii and only if there exists a reaction 1^ — >■ 
Si + Su for some Su £ 5^- We have bij 7^ if and only if dj^i 7^ or 6^^^^ 7^ for some 
k and t, and hence either there is a reaction Sj — ;■ Si or there exist Su,St G S^, so that 
Sj+ St ultimately reacts to Si+Su via y^. If Gy^ is strongly connected then by Proposition 
|3.1[ ii) the existence of these reactions is a sufficient condition. It follows, after reordering 
of the species in Sa , that S is a block diagonal matrix with blocks given by the species in 
each connected component of Gsa,ya- Indeed, if Si,Sj are in different components, then 

It follows from the lemma that a necessary condition for bij 7^ is that Si can be "pro- 
duced" from Sj . We restrict the study to the case where Gs^ ,ya is connected and note that 
the results apply to every connected component individually. However, the propositions to 
be derived below are stated in full generality, that is, without the assumption that Gs^ya 
is connected. 



Using (3.8), the column sums of B are zero. Thus, B is the Laplacian of a labeled 
directed graph Gs^ with node set Sa and an edge from Sj to Si whenever bij 7^ 0, z 7^ j. 
Note that lij £ M(Con)[5^] is S-positive. 

Since Gs^y^ is connected, then so is Gs^- In general, two species Si,Sj belong to the 
same connected component of Gs^y^ if and only if they belong to the same connected 
component of Gs^ . We will use this fact repeatedly in what follows. 

By the Matrix- Tree theorem, the principal minors B^j) oi B = C{Gs^) are 

TGOiSj) 

Thus, B has rank Na — 1 if and only if there exists at least one spanning tree in Gs^ rooted 
at some Sj with j G {!,••• ,Na}. For a general PTM system with a selected cut Sa, we 
obtain the following proposition: 

Proposition 3.10. The non-interacting graphs provide all conservation laws involving 
only the substrates Sa,i if and only if Gs^ ; has at least one rooted spanning tree for all I. 

Proof. The non-interacting graphs provide all conservation laws involving only Sa^i if and 
only if all conservation laws are multiples of Yls i^s ^i ~^ Sy ey ^k = 0, which is the 
case if and only if the rank of Bi is Na^i — 1- As stated above this is equivalent to the 
existence of a rooted spanning tree in Gs • D 
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Remark. In particular, the lemma holds if Gs^ is strongly connected. If Gy^ is strongly 
connected, then to check that Gs^ is strongly connected we do not need to calculate the 
labels of Gs^ ■ Whether there is an edge or not between two nodes follows from the set of 



reactions, cf. Lemma 3.9 



For simplicity we assume that there exists a spanning tree rooted at Si. Then, the 
variables 5*2, . . . , Snc, can be solved in the coefficient field M(ConU5^U{5i}). In particular, 
using Cramer's rule and the Matrix- Tree theorem, we obtain 
(3.11) 



S. 



' "'" = -if^Si = r {S^^)Si, where ^ " ^^ 



"(1,1) '^(Sg)- •^^""'~^' ••"'"" icx,(5S) =E..e(5,)^W 

and j = 2, . . . .,Na. It follows that o"(5^) is S-positive and (Tj{S%) is either a zero or S- 
positive element of M(Con)[5^]. If the graph Gg^ is strongly connected, then (Tj(5^) / 
for all j and any choice of Sj could be used instead of 5*1. Further: 

Proposition 3.12. A connected component Gg^i of the graph Gs^ is strongly connected 
if and only if aj{S^) is a non-zero rational function in ]R(ConU5^) for all Sj € Sa,i- 

The results shown above provide a proof of the following lemma. 

Lemma 3.13. If a substrate St ^ S^ is a variable in the rational function r^{S^) for some 
Sj £ Soi, then there is Si G Sa and Su G 5^, such that Si + St ultimately reacts to Sj + Su 
via ya ■ 



After substitution of the value of 5*^ given in (3.11) into Y^ (3.2) we obtain 

(3.14) Yk = rl{S^)Si, 

where r^ is either zero or an S-positive rational function in S^ with coefficients in M(Con). 
If Gy^ is strongly connected then this function is non-zero. 

Conservation laws. The sum of the species concentrations in Gs^y^ is conserved. If 
the total amount S*! = S*! + • • • + Sn^ + ^i + • • • + Yp^ is given, we obtain 

S, = {1+ rf(5^) + • • • + 4J5^) + rliS^) + ■■■ + r^^iS^) )5i, 

where the coefficient of Si is an S-positive element of M(ConU5^) and thus, 

with ff an S-positive rational function in S^ with coefficients in ]R(ConU{5i}). 

Further, if 5i > then S*! 7^ at steady state and 5i > for non-negative values of the 



substrates in S^. This remark and Proposition 3.12 imply: 

Proposition 3.15. A connected component Gg^i of the graph Gs^ is strongly connected 
if and only if any steady state solution satisfies Sj 7^ for all Sj E Sa^i, and any total 
amounts Si > 0. 

By substitution of 5i by ff, we obtain 

(3.16) Y, = rliSl) := r,^ (5^)rf (5^), S, = rf (5^) := rf (5S)ff (5^^ 



■^aJ 



with r\ ,rf either zero or S-positive rational functions in 5^ with coefficients in M(ConU{S'i}). 
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Proposition 3.17. Assume that for each I = l,...,nQ,, there exists a spanning tree of 



Gs^i,y^i rooted at some species Si^. Then, equations (3.5) are satisfied if and only if 

where r^ is zero or an S-positive rational function in S^ with coefficients in M(Con). Fur- 
ther, the conservation law Si = J2s es '^i ~^ "l^Y-ey ^k is fulfilled if and only if 

(3.18) S,,=ff^{S'J, 

where rf is an S-positive rational function in S^ with coefficients in M(ConU{S';}). 



In Example (2.1), the graph Gsi has two connected components: 5*5, which does not 

^2,1 

ahow further ehminations, and 5i <_ > S2 , which is strongly connected. Selecting Si as 



f-i,: 



the non-eliminated species we obtain 

Q ^2,1 c, ^1.2 c- V 
J2 — ~ — '-'1 — 79 9 — TT'-'l; -'I - 



^2,1 g ^ ^1.2 

^1,2 ^1,4^2,3'S'3 



79 9— '-'1 — 79 "Jl) J-2 



,3 



6? 



The total amount equations Si 

'Sl = s^{l + ^il^s^), 



'1,4^1,2 

55 + Fs and ^2 = Si + S2 + Yi+ Y2 give: 

di2 i 1 &?4 + C2,l , h\ 



S2 



^J.US3 



+ 



+ 1 + 



Cl,2 



dl,2 



^1,2 c- 
"1,4 



Si. 



Letff(53,S4) = 52 



fii-^Sa 



+ 



"U 



-^^ + 1 + 5^ 



-1 



thus: 



(3.19) 



Si 



-j^rf{S3,S4) 
"1,2 



^2 



rf(g3,g4) 

/^2,3'5'3 



(&f,4 + C2,l) ^g 

Fi = — ^ rf(53,54) 

Cl,2 



y2 = rf(S3,S4), 



5. 



X3 



Si 



1 + /u| 554 ■ 

/^4,5'S'l'S'4 
1 + /XI554 ■ 



Thus, all species are given as S-positive rational functions of ^3,54 in the coefficient field 

M(ConU{5i,52}). 

3.4. Steady state equations. To summarize, at steady state the intermediate complexes 
y can be expressed as rational functions of the substrates S and therefore eliminated. Fur- 
ther, provided a cut Sa exists, the variables Sa can be expressed as functions of 5^ = S\Sa 
and therefore also eliminated. For the latter statement, we make use of the conservation 
laws (with given total amounts) for the species in Sa determined by the connected compo- 
nents oiGs^,y^. 



Specifically, consider the steady state equations (3.5 ) for 5^. Substituting the expressions 



in (3.14) and (3.11) for 3^ and Sa provides the steady states equations in terms of Sa and 
the selected variables S^ (one for each conencted component of Gs^). Using (3.16), the 



steady states equations are given in terms of 5^ only. Let 5^ = {Sn^+i, . . . , Sn} and 
let ^u{<Sa) = be the equation obtained from Su = after elimination of y and Sa and 
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removal of denominators. The denominators can be cfiosen to be S-positive and we can 
multiply tfie expressions by the denominators without changing the positive solutions. 

Assume that the graph Gs^y^ has ria connected components, and recall that each 
of them gives rise to only one conservation law (Proposition 3.10). Extend the set of 
conservation laws to a maximal set of dim(r-'-) laws. 

Theorem 3.20. Consider a PTM system for which there exists a cut Sa- Further, assume 
that each connected component of Gs^ admits a rooted spanning tree. If total amounts 
Si are given for the Ua connected components of Gs^ya ^^'^ ^^^ dim(r ) — Ua additional 
conservation laws, then the non-negative steady states of the system with positive values 
for all substrates in S^ are in one-to-one correspondence with the positive solutions to 



$„(5S) = 0, Si = ^i{Sl 



for u = Na + 1, . . . ,N and I 



rir 



+ l,...,dim(r^ 



Proof. We have shown that any non-negative steady state solution with positive values 
for all substrates in S^ must satisfy these equations. For the reverse, consider a positive 
solution s = (sat^+i, . . . , Sat) to the equations <I>u(5^) = and Si = Lpi{S%). For i = 
1, . . . ,Na, define Si through equation (3.18) and yk, k = 1, . . . , P, through equation (3.2). 



For positive rate constants and positive total amounts, Si,yk are non-negative (because of 
the S-positivity of the rational functions defining them) . By construction these definitions 
automatically ensure that the conservation laws with total amounts Si, I = 1, . . . , Na, are 



satisfied (see Proposition 3.17). 



By Proposition 3.1 the values yi,...,yp satisfy (2.2) for all k and hence the steady 



state equations of the intermediate complexes are satisfied. By Proposition 3.17 the values 
si, . . . , SNa satisfy (3.5). Since the latter is just (2.3) after substitution of (3.2), we see that 



(2.3) holds as well. Since <I>u(5^) = is the steady state equation S'„ = after substitution 
of (3.2) and (3.18), this equation is also satisfied and the same reasoning applies to the 
equation Si = (/?;(5^), I > Ua- Thus, Si = Si and Yk = yk is a solution to the steady state 



equations and satisfy the conservation laws corresponding to the total amounts Si. 



D 



This theorem together with Proposition |3.2[ iii) and Proposition 3.15 gives the following 
corollary. 

Corollary 3.21. Assume that Gy is strongly connected and that for all S £ S there exists 
a cut Sa such that S £ Sa and Gs^ is strongly connected. Then, Si = or Y^ = is not 



a steady state solution for any i,k. With the notation of Theorem 3.20, the non-negative 
steady states of the system are in one-to-one correspondence with the non-negative solutions 
to 



for u = Na + 1, 



$„(5S) = 0, Si = ^i{Sl) 
N and / = no + 1, . . . ,dim(r ). 



In Example (2.1), dim(r ) — n^ = 1 and only one conservation law is missing, Sa = 
Sa + 5*4 + li + 12 + ^3- The elimination procedure leads to the steady state equations 
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consisting of 5*3 = {^3) and Ss {(ps): 

^3,5^4,5 •S'i5'4 



MS3,S4) = -bt,?US3,Si) + 



1 + ^ii.S^ 



S3 = MS3, Si) = S3 + S, + -^ ■ ff (53, 54) + "^''^ 



ci,2 1 + ^J'%^^s. 



4 



Since tfie conditions of Corollary 3.21 are fulfilled, any non-negative solution of this reduced 



system provides a non-negative steady state of the PTM system. The steady states of the 



other species, 5*1, S2, 5*5, Yi, Y2, Y3, are found from (3.19). In this specific example, the first 
equation is easily transformed into a linear equation in ^3,^4, and hence either 5*3 or ^4 
can be eliminated as well, providing a polynomial equation in the remaining variables. In 
this case, S-positivity is not guaranteed. 

In the example we intentionally selected Sa to have the highest possible number of 



elements, since all these variables are subsequently eliminated. In Example (2.10), the 
cut Sa = {E, S3, Si} allows us to eliminate three substrates and reduce the steady state 
equations to a system of three equations in three variables. 



In some systems (see e.g. Section 4.2) there two different cuts Sa, S'a might exist, such 
that the union is not a cut, but still all variables in 5q, U5^ can be eliminated. Thus, more 
species might be eliminated if different cuts are considered. 

4. Examples 

4.1. TG frame'work. In [TTj, the authors provide a linear elimination procedure for the 
special case in which the set of substrates is partitioned into two distinct sets. In their 
context, a PTM system (here called TG system) consists of three non-empty and disjoint 
sets of species called enzymes, substrates, and intermediate complexes: 

Enz = {El,..., El}, Snh = {Si, ..., Sn}, Int = {Yi, . . . ,yp}, 

and a set of reactions Ret = Ra U RbU Re with 

Ra = {Ei + Sj ^Yk\{i, j, k)ela} Rc = {Y,^Yj\{i,j)e Ic} 

Rb = {Yk^E, + Sj\{i,j,k)€h} 

for /„, lb C {1, . . . , L} X {1, . . . , iV} X {1, . . . , P} and Ic C {1, . . . , P}2, such that (i) All 
chemical species are involved in at least one reaction; (ii) For every intermediate complex 
Yk there is at most one enzyme i5^^(fe), such that {rj{k),j, k) £ RaU Rb for some j; (iii) If 

two intermediate complexes Yk,Yy are 1-linked, then £'^(fc) = E^^i^^y Further, the graph Gy 
and each connected component of the graph Gsub are required to be strongly connected. 
In particular, the assumption that Gy is strongly connected implies that any 1^ ultimately 
reacts to Si + Sj for some i,j. This is our Assumption (ii) of a PTM system. 

Essentially, they consider post-translational modification systems in which the enzymes 
are not allowed to be modified. Let S = SubUEnz, Sa = Sub, and 5^ = Enz. Properties 
(i)-(iii) imply that Sa is a cut. Note that ya = ya — 3^- Thus the framework developed 
here is an extension of the framework developed in |17] . 
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By assumption (iii) the graph Gsgy has L connected components that provide L con- 
servation laws for the enzymes: Ei = -Ej + ^^.i (^wj Ifcj for i = 1, . . . , L. With the notation 



of Lemma 2.13, A''^ = n'^ = L, P^ = 0, so that A^^ + P^ — n^ = and thus a set of 
independent conservation laws of a TG system can be derived from the non-interacting 
graphs of Gs,y- Further, the form of Ra and i?^ ensures that any non-interacting graph 
contains species either from Enz or Sub, but not both. Thus, all conservation laws are 
associated with a connected component either of GEnz,y or Gsub,y- 

It follows that if all intermediate complexes ultimately dissociate into an enzyme and 
a substrate, and each connected component of Gg^ admits a rooted spanning tree, then 
elimination of the variables in 5q, U 3^ reduces the steady state equations to L equations 
derived from the total amount of enzymes. 

4.2. Signaling cascades. Our setting is well-suited to study elimination of variables in 
signaling pathways. Signaling pathways form a special type of PTM systems and an ex- 
tension of TG systems to include some substrates that also act as enzymes. 

Definition 4.1. A signaling cascade is a collection of TG systems R^, . . . , i?", with corre- 
sponding sets of species 

Enz^ = {E{,...,El}, SnW = {Si,..., S'^J, T = {yI, . . . ,yI,J 

and sets of reactions Ret* = i?^Ui?^Ui?*, for i = 1, . . . ,n, satisfying the following conditions: 

(i) (Enz* U Sub* U y^) n (Enz^' U SuW U y^) = {Ei^^} = {Sj^j iij = i + l and it is empty 

otherwise, 
(ii) For all i, each connected component of the graph Gg^j^i admits a spanning tree rooted 

at S}^_. 
(iii) All intermediate complexes ultimately dissociate into two substrates. 



Condition (i) implies that a signaling cascade consists of independent TG systems 
"joined" by only one substrate acting as an enzyme in the layer below. This descrip- 
tion fits signaling pathways in which the signal is transmitted downstream. Condition (ii) 
ensures that the intermediate complexes can be eliminated. 

Let A^ = A^i H \- Nn, L = Li ^ \- Ln and 5 = Uj Enz* U Sub*. For each i, consider 

the subset Sub* C S. The associated set of intermediate complexes is 3^sub» = 3^* U {Yk G 
y^~^^\ rj{k) = S]^.}, and Sub* is closed (TG systems do not incorporate reactions Su — >• 5*^ 
among substrates or enzymes). By definition, substrates in Sub* do not interact and thus 
Sub* is a cut. 

For simplicity, we assume that the graph Gg^^^i is connected for each i. By Proposition 



3.17, elimination of the variables in Sub* provides the steady state relation 

S;*=r*(Enz*)5]v,, 5*GSub*\{%}. 



By Lemma 3.13, r* depends on the species in Enz* only: if S"^ -|- St ultimately reacts to 



S* -|- Sr for some species 5^ in Sub* and Sr ^ S \ Sub* via y , then since 5*- ^ S^y^, 
St = Sr = Ej^ for some Ej^ £ Enz*. Further, if Yk £ 3^subM ^^ let Yfc = ^^(Enz*)^^. be the 
corresponding rational function. 
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Conservation laws. Since Gg^j^j is connected and admits a rooted spanning tree, the 
sum of tfie species in tlie grapli Gg^j^i y , provides the only conservation law among the 

species in Sub* Uj^g^j^i. Thus, for each i, let a total amount Si be given. We have at steady 
state 

(4.2) S,= Yl ^i(Enz*)5k+ E ^r(Enz^)^k- 

For i = n, S'^ ^ Enz", and so S"^ is expressed as a rational function in Enz". 

Thus, if we let Enz = IJ^Enz*, we have that the species in 5 \ Enz are given as 
rational functions in Enz with coefficients in M(Con). Condition (iii) implies that for 
E £ Enz* VJS'^y,}, {E} is a cut with associated (connected) graph GEy^- Thus, if the total 
amount E is provided, the steady states must fulfill the equality 

(4.3) E = E+ Y. ^k = E+ Y ^r(Enz*)%. 

We conclude that the non-negative steady states of a signaling cascade are solutions to L 
equations in Enz with coefficients in M(Con), provided that total amounts for Enz are given; 
that is, Si,...,^^^! for the enzymes S}^^, (4.2) and E^^ for fi^ S Enz\{S^^, . . . , S']^~_ }, 



(14J. 

Note that the number of conservation laws obtained in this way is m = ^^ Lj + 1 
(remember Sn)- Let e = 1 if n is even and otherwise, and let e^ = 1 — e. The cuts provide 
all conservation laws: The graph associated to the cut 



5„ = y Sub' U y Enz* 



i odd 



has Ua = e + Yli odd -^« connected components and thus, n'^ = e'^ + ^^ ^^^^ Li. We have 

Na = Ei odd L^ + E^ even(^^ " 1) + f, and N^ = Z^ even ^« + E^ odd(^^^ " 1) + ^'- Farther, 

y^ = y^ so that P^ = 0. 



Let 7 = N^ ~ "-a = Yli odd(-^« ~ ^)- -^y Lemma 2.13, if there are 7 independent terms 
in S^ n r, then all conservation laws come from non-interacting graphs. By hypothesis, 
for i even, the graph Gg^j^i has a spanning tree rooted at some node Sj. This means that 
for every Su / Sj in Sub*, there exists a directed path Su -^ 'S'fci ^ ■ ■ ■ -^ Sk^ ^ Sj . By 



the conditions of a TG system and Lemma 3.9, an edge S^^ — >• S^^ implies that there is 



either a reaction Sk^ — t- S^^ , or E + S^^ ultimately reacts to i? -|- 5^^ via y. In either case, 
we see that Su — Sj £ S^CiT for all Su / Sj in Sub*, implying that there are indeed 7 
independent vectors in 5^ n L. 

4.3. Biological examples. 

MAPK signaling cascade. We consider the first two layers of the MAPK cascade: a 
two-layer cascade with one-site modification in the first layer and two-site modifications in 
the second layer. In the latter, dephosphorylation is considered sequential but this is not 
the case for phosphorylation |15j. 
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The reactions of the system in the first layer are 

E + Sl^^Yl^E + Si Fi + Si ^± Yi ^ Fi + S^ 

accounting for phosphorylation and dephosphorylation, respectively, via a Michaelis-Menten 
mechanism. In the second layer we have the phosphorylation reactions 

'-'1 + '-'0,0 -« — ^1 ^ '-'1 ^ '-'1,0 '-'1 ^ '-'0,0 -« — ^2 '-'1 ^ '-'0,1 

c-l I q2 — s. -1^2 . cl I C2 cl I c<2 — ^ -1^2 . cl I c2 

'-'1 + '-'1,0 -« — J^3 ^ '-'1 ^ '-'1,1 '-'1 ^ '-'0,1 -^ — -^4 ^ '-'1 ^ '-'1,1 

Dephosphorylation proceeds sequentially in the following way: 

F2 + Si, ^ Yi -^F, + Si, F2 + Si, -r Yi -^F, + Si, 

The sets of enzymes are Enz = {E,Fi}, Enz = {Sl,F2}. The sets of substrates are 
Sub^ = {5'o,5|}, Sub^ = {S,,, SiQ,SQi,Sfi}. The sets of intermediate complexes are 
Int^ = {Yl,Yi}, Int^ = {Yi,Yi,Yi,Y^,Yi',Yi}. We have Enz^ n Sub^ = {Sl}, so that 
the modified substrate in the first layer is a kinase of the next layer. The superindex 
denotes the layer, while the subindex denotes phosphorylation state (the presence of the 
phosphate group is represented by 1). 

The components of the graph Gy are each of the intermediate complexes and are thus 
strongly connected. The graphs G'3y^ji and Gg^{j2 are 



•^0 < ^ ^1 =r 1. 



c2 >=' \ C2 

'-'0,0 ^ _ '-'1,1 



=^02" 
'-'0,1 

which are also strongly connected. The conservations laws (all derived from non-interacting 
graphs) are 

E = E + Y,^ Si = sl + Sl + y/ + Yi + Yi + Yi + Yi + Yi 

Fi = Fi + Yi S2 = Slo + Slo + Sl, + Sl, + Yi + Yi + Yi + Yi + Yi + Yi 

F2 = F2 + Yi + Yi 

Therefore, if total amounts are provided, then the steady states of the two-layer cascade 
are found as solutions to a system of four polynomial equations in four variables, namely 
-E'i-^i)-^2,'S'i- 

Receptor protein-tyrosine kinase. Receptor protein-tyrosine kinases (RPTK) are 
cell surface receptors linked to enzymes that phosphorylate their substrate proteins in tyro- 
sine residues. The common mechanism for their activation is autophosphorylation following 
ligand-induced dimerization [21 §15]. The phosphorylated receptor serves as binding site 
to downstream signaling molecules, such as SH2 domain containing proteins. Further, the 
receptor can be dephosphorylated by several protein tyrosine phosphatases (FTP) |16j . 

A simple model describing the phosphorylation state of an RFTK is: 

2Ro^^Yi^2Ri S + Ri^^Y2 F + Ri:^Y3^F + Rq 
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where Rq,Ri stands for the unphosphorylated and phosphorylated RPTK respectively, S 
is a protein binding Ri, and F is a PTP. 

We have S = {Ro,Ri,S,F} and y = {Yi,Y2,Y3}. Note that Sq = {Ro,Ri} are the 
self-interacting substrates and thus cannot be part of a cut. First of all, the intermediate 
complexes Yfc can be eliminated in terms of S. The graph Gsy is 



Rn 



F — Ys 



Y1-R1-Y2 — S 



The non-interacting graphs provide two conservation laws: F = F + Y3, and S = S + Y2, 
associated to the cut Sa = {F, S}. Thus, the substrates F, S can be eliminated, in fact from 
the conservation laws. We conclude that at steady state all species are described as rational 
functions of Ro,Ri and the non-negative steady states are in one-to-one correspondence 
with the non-negative solutions to the equations corresponding to Rq and the remaining 
conservation law R = Rq + Ri + 2Yi -1-1^-1-13. 
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